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Abstract 

A kink-based path integral method, previously applied to atomic systems, is modified and used to 
study molecular systems. The method allows the simultaneous evolution of atomic and electronic 
degrees of freedom. Results for CH4, NH3, and H2O demonstrate this method to be accurate for 
both geometries and energies. Comparison with DFT and MP2 level calculations show the path 
integral approach to produce energies in close agreement with MP2 energies and geometries in 
close agreement with both DFT and MP2 results. 
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I. INTRODUCTION 



The development of simulation methods that are capable of treating electronic degrees 
of freedom at finite temperatures is necessary for the study of a variety of important 
systems including those with multiple isomers with similar energies (such as metal clus- 
ters) and with dynamic bond breaking/forming processes. A fundamental difficulty in 
using ab initio quantum approaches to study systems at finite temperatures is the need 
for most algorithms to solve a quantum problem (to find, for example, the ab initio 
forces) at each geometric configuration. Thus the CPU requirement per time or Monte 
Carlo step often prevents a simulation. Feynman's path integral formulation of quantum 
mechanics offers the possibility of simultaneously treating geometric and electronic de- 
grees of freedom without the restriction of solving a quantum problem for fixed atomic 
positions. In addition, temperature and electron-electron correlation can be included and 
make this approach very tempting as a starting point for ab initio simulations. An un- 
fortunate aspect of the path integral approach is the so-called "sign problem" which can 
make the standard deviations of estimate d q uantities (such as energy) too large for prac- 
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problem occurs because the quantum density matrix is not positive definite and results 
in averages being determined from sums of large numbers with different signs. The Car- 
Parinello implementation of density functional theory [2J] is motivated by needs similar to 
the ones described above and treats electronic and geometric degrees of freedom on a similar 
footing and allows both types of degrees of freedom to propagate simultaneously during a 
calculation. A limitation of this approach, though, is the need for the electronic degrees of 
freedom, as described using single particle orbitals, to be very close to the lowest energy set 
of orbitals, forcing the use of small time steps in a molecular dynamics simulation. 

We have recently introduced 0, Q a "kink-based" path integral approach and have 
demonstrated that it can be used to overcome the "sign problem" in atomic systems. In 
the present work a formalism appropriate for molecular systems is developed. To construct 
a practical approach, an approximation to the exact path integral approach is made; the 
approximation is based on the results of our previous work. The method treats the electronic 
structure as consisting of ground and excited single determinant states built from atom- 
based orbitals. Simulations include moves that perform unitary transformations of the 
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single particle orbitals, additions and deletions from a list of excited states used to evaluate 
the energy, and moves of atoms. Using this procedure, electronic and geometric degrees of 
freedom are treated simultaneously. The method is used with success to calculate the average 
energies and geometries of CH 4 , NH 3 , and H 2 at finite temperatures. We define success as 
(a) overcoming the sign problem, (b) not requiring low energy orbitals, (c) average molecular 
geometries in agreement with previous ab initio calculations, and (d) average energies that 
compare favorably with previous ab initio calculations. 



II. KINK-BASED PATH INTEGRAL FORMULATION 

Our previous workjiol. I^J started with the path integral expression for the canonical 
partition function, evaluated with fixed geometries and in a discrete N-electron basis set { | j > 

}: 



Q({\j >}) = Tr{exp(-(3H)} = £ < j\ exp(-/?F)|j > 



(1) 



Making the Trotter approximation and discretizing the path into P segments, we find 

Q(P,{\j >}) = (TIE) < Mexp(-I3H/P)\j i+1 > (2) 

\i=l ji / 

This can be interpreted as a path in the space of states that starts and ends with state 
\ji >. If \ji >= >, we have a diagonal matrix element, otherwise we have an off- 

diagonal matrix element. Any place that an off-diagonal element appears is called a "kink" 
and it is clear that the paths can be classified into those paths with zero kinks, 2 kinks, 3 



kinks, etc. We demonstrated 
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(3) 



(4) 
(5) 



(3 = 1/fcgT, P is the discretizing variable in the path integral formulation, n is the number 
of kinks, Sj is the number of times a particular state appears in the sum, and m is the 
number of distinct states that appear (YlT=i s j — n). As written, Eqn. 01 is amenable to a 
Monte Carlo simulation. However, since the matrix elements can be negative, the usual sign 



problem will occur if the states are not well chosen. In our previous work [19, |20j], the initial 
N-electron states \j > were chosen to be simple, anti-symmetrized products of 1-electron 
orbitals. The N-electron states were improved by periodically diagonalizing the Hamiltonian 
in the space of those states that occurred during the simulation. The result was that the final 
states were linear combinations of the initial orbitals (essentially the "natural spin orbitals" 
for the system), the only paths that appeared at the end of the simulation contained 0, 2, or 
3 kinks, and the sign problem was reduced to insignificance. Further, the dimensions of the 
density matrix were small enough that the matrix could be kept in memory and transformed 
whenever a diagonalization took place; this did not require any matrix elements involving 
the initial orbitals, which significantly reduced the computational effort. 

In the case of geometric degrees of freedom, all matrix elements are referenced to the 
initial orbitals and the adaptive scheme used for atomic systems must be modified for com- 
putational efficacy Our previous work showed that once a good guess for the ground state 
was obtained the vast majority of paths contained or 2 kinks. Using an approximate 
Hartree-Fock solution as a guess for the ground state, an approximate infinite order summa- 
tion of kinks from the ground state is developed (leaving for future work the straightforward 
extension to the case of degenerate or nearly degenerate ground states), in which we assume 
the most important paths contain many instances of the ground state. 

First consider a Hartree-Fock-like approximation to Q{P). In the Born-Oppenheimer 
approximation, Q is given by 

Q = Tr {exp (-(3H)} = f dR N £ < j\ exp {-(3H(R N ) \j > (6) 

J j 

where we assume N nuclei, N e electrons, and that {\j >} is a set of N e -electron orbitals. Each 
N e -electron orbital is expressed as an anti-symmetrized product of 1-electron orthonormal- 
ized spin-orbitals \j >= A^^a^, <pj 2 aj 2 , (pj N aj N ) and typically these 1-electron spatial 
orbitals are expressed in terms of atom-centered orbitals {\x% >} (themselves often sums or 
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"contractions" of gaussian orbitals) 

I0j > = J2 c v\ Xi > ( ? ) 

i 

and the Hamiltonian matrix elements are expressed in terms of {\xi >}■ For a given geome- 
try, the Hartree-Fock orbitals will be a unitary transformation of any arbitrary starting set 
of orbitals {|0° >} 

\^ F > = > ( 8 ) 

3 

\j HF > = A{^ F a h ^f F a j2) ...) (9) 

Symbolically, we will write this as \ j HF >= \U HF j° > and an arbitrary unitary transforma- 
tion as \j >= \U j° >. Since the trace is invariant with respect to unitary transformations, 

Q = [dR N Y,< f\ exp(-pH(R N )\f > (10) 

J 3 

oc /"dR w ^^<f/j |exp(-^(R JV )|f/j°> (11) 

J U jo 

where the proportionality constant depends on the number of possible unitary transforms. 
Treating the allowed unitary transforms as rotations, the proportionality constant will then 
be a proportional to a power of tt and the sum over U becomes an integral over rotational 
angles. As written, it is clear that a possible algorithm is to view the unitary transforms 
as just another degree of freedom to be sampled in a Monte Carlo or molecular dynamics 
simulation. 

This expression is exact within the use of a finite set of states and any related basis 
set superposition errors. To make progress, we will make approximations to the density 
matrix elements < j|exp [— /^(R^)] \k >. The most obvious one is a Hartree-Fock-like 
approximation 

< j\ exp [-f3H(R N )} \k > w exp [-^-(R^)] 5 jjk (12) 
Eqn. ^2 then becomes 

Q HF cc fdR N ^2^2exp[-/3H(R N ) Vj0tUj o] (13) 

J u jo 

The delta function in Eqn. ^] in essence generates paths with kinks. Thus, paths with 
kinks will correspond to a Hartree-Fock approximation and no electron correlation will be 



included in a calculation using Eqn. ED This expression is an approximation to Eqn. ^2 
in two respects; the obvious one being the approximation to the density matrix and a less 
apparent one that the sum is not independent of U. In terms of a simulation in which 
U is sampled, there is an entropy associated with the different U's which means that not 
just U HF will appear, but other U's, which may result in average energies that are higher 
than the Hartree-Fock energy. Of course, sophisticated simulation methods can be used to 
minimize the effects of this entropy. 

To go beyond the Hartree-Fock approximation and include correlation, consider the dis- 
cretized version for Q, Eqn. Eland write it in the Born-Oppenheimer approximation as 



Q(P) = |dR^jj>f(R*)+ 
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As a first step, we assume that the most important paths with at least 2 kinks will consist of 
alternating ground and excited states. That is, half of the states will be the ground state and 
the other half will be excited states. Assuming as previously stated that the lowest energy 
state is non-degenerate, we can write (where now the summation variable n denotes the 
number of times the lowest energy state appears in a path and we suppress the dependence 
of Xj and £^ on H, N for notational convenience) Eqn. as 
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where the binomial factor accounts for the number of different ways to make the different 
excited states appear in the path and the factor of 2 appears because the first state in the 
path can be either the ground or an excited state. Empirically, we have found the the 
most important term in the sum over m is the m = term. This can understood from the 
following; the ratio of the m = 1 term to the m = is 

(n — l)xonFi 



(P-n + l)r 

where 



(19) 



dXQ \%0 Xj 

Now for small values of (3/P = e we have 



x = < 0|exp(-eiT)|0 >^ 1 - eE (21) 

t oj = < 0| exp(-e#|0 >^ -eH 0j (22) 
H 2 

^ ~ 6 E F 0J 77 = - eAi? MP2 (23) 

|r./r.| « ^ (24) 



(n — l)x nTi n(n — 1) 



(25) 



(P-n+ l)r /3A£ 

where AEmp2 is the MP2 correction to the Hartree-Fock energy and AE represents a typical 
difference in energy between the lowest energy state and one of the excited states appearing 
in the sum for r . We typically expect (3AE » 1 and hence the m = term to be the 
most important in Eqn. ^] as long as the sum on n is quickly converging. 

Returning to Eqn. we evaluate just the m = to find the interesting result 

Qf = f dR N { x? + £ ( ) (26) 



n=l 



dR iV {(a;o + ro) p } (27) 

where we assume that the sum on n converges sufficiently quickly so that the sum can be 
extended from P/2 to P; this accuracy of this assumption was checked and confirmed in the 
calculations performed in this paper. Further insight can be gained when e is very small. In 
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this case 

Qf « J dR N {(1 - e (£ + A£ M P 2 )) P } (28) 
< £ >= -d\nQ/d(3 « < £ + AEmp2 >rn=< E MP2 > nN (29) 

where the subscripts on the last averages indicate averaging over the geometric degrees of 
freedom. So paths of alternating ground and excited states, when only the m = term is 
included, should be expected to give an MP2-level result. Eqn. EH will be accurate when a 
very good guess to the lowest energy state corresponding to the Hartree-Fock solution exists 
for a given geometry. However, this will not always be the case during a simulation and 
Eqn. will be refined by including terms beyond the m = term and more complicated 
kink patterns in which there is more than a single excited state between occurrences of the 
ground state. 

We first consider the m > term in Eqn. ^] and consider more complicated kink patterns 
later. Note that m > will include Ti, the first derivative of T , and higher order derivatives 
of IV We expect, and have verified in the systems we have studied, that the major correction 
to Eqn. EH! contains terms with only IV The derivation of an expression that includes 
derivatives up to second order, is possible, but not included in herein. Including only those 
terms in Eqn. E3of l ess than second order in derivatives of Tq, we find 

2 P ^ (n - 1\ [P - l)\x?- n+m n!r™-"T7 
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So we can then find 

P/2 



Q? = J dR N | 



k=l 



»+E(r>o M a h/ dRw <- + t^< 32 > 



To include paths with more than one excited state between each occurrence of the lowest 
energy state, the above approach is generalized. We will refer to a portion of the path 
between 2 occurrences of the ground state as an "excursion". An excursion will contain one 
or more excited states. The development of an expression for the "lowest energy dominated" 
(LED) set of paths begins by defining a "weight" associated with any particular excursion 
j to be 

= fo, -*.)(*„- fo, -x.) <33) 
where the excursion j is defined to include the excited states a,b, • • • ,z. Then 




n=l v ' u nino u i 



where 2n + An({n.j}) is the number of kinks for a particular set of excursions. If the states 
are well chosen, we expect the contributions from excursions with greater than one excited 
state per excursion to be much less than the contributions from the one excited state per 
excursion set of paths. Thus, as an approximation, we write 

1 1 



2n + An({n i }) ' ' 2n 



(35) 



and we find 



Q LED ~ J dR- 4 + E S^" lr ° (36) 



where 



j X 3 ( X - Xj) (X - X k ) 



Defining two matrices, W and M 



(Woh = t0tt ° J = (38) 

y/(x - Xi)(x - Xj) 

(M ) y = \ — (1 - 5 id ) (39) 

a/^o - xj(x - Xj) 



obtains 



r = Tr {W + W ■ M + W ■ M ■ M + ■ ■ ■ } = Tr {W ■ (I - Mo)' 1 } (40) 

This sums to infinite order all possible types of excursions from the lowest energy state, 
with the proviso that the contributions from excursions with different numbers of states is 
a rapidly decreasing function of the number of states involved in the excursion. With 

r\ = -^-r = Tr {W 1 ■ (I - Mo)' 1 + W ■ (I - Mo)' 1 ■ M 1 ■ (I - Mo)" 1 } (41) 
dxo 

we immediately find 

Q LED = /dR^{(x + ^) P }oc|dR^^{(x ( f /) + r I^r} (42) 

Eqn. 1321 is the principle result of this work and represents an expression that can be used 
to simulate a molecular system. There are two computational challenges to using this 
equation. The first is that the sum over excited states includes all states and can become a 
severe bottleneck in a calculation. This issue can be addressed using a limited set of excited 
states such as 1- and 2-electron excited states (an approximation made in this work) and 
by making a slight modification to Eqn. 021 that will enable excited states to be sampled 
during the Monte Carlo process. Qled is a function of the ground and excited states, 
whose total number can be very large. However, it is expected that only a subset of the 
states will contribute significantly to the partition function and thus we wish to develop 
a Monte Carlo procedure that will limit the number of excited states used to those with 
significant contributions to the partition function, as judged by a Monte Carlo simulation. 
To do this, first label the excited states in order of decreasing magnitude of [^/(xo ~ 
(approximately an excited state's contribution to the MP2 energy). Next the excited states 
are divided into groups of N g states, with group 1 corresponding to excited states 1 . . . N g , 
group 2 to excited states N g + 1 . . . 2N g , etc. If there is a total of M g such groups and 
QhEoiP'g) denotes the result obtained using only the first n g groups of excited states, Qled 
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becomes 



QlED — QLED^Mg) 

= Qled(M 9 ) - Q LED (M g - 1) + Q LED (M g - 1) 

= [Qled(M 9 ) - Q LED (M g - 1)] + [Q LED (M g - 1) - Q LED (M g - 2)] + 

■■■ + [Qled{1) -Qhf]+Qhf 
= AQ LED (M g ) + AQ LED (M g -!) + ■•• + AQ LED (l) + AQ LED (0) 



In this notation AQ EEE (0) = Qhf and j can be sampled during the Monte Carlo procedure. 
If the states are reasonably ordered, the sum over j should converge for a relatively small 
value of j and the matrices involved in evaluating AQl E dU) will be manageable in size. 
The other time consuming part of any calculation will be the inversion of the matrix I — M . 
Fortunately, this is amenable to parallel computation using standard linear algebra packages 
which will aid in the implementation of the method. 

III. MONTE CARLO SAMPLING PROCEDURE 
A. Rotations of single particle states 

Sampling unitary transformations U was accomplished in the following way. Two orbitals 
<f>j and <pi were randomly chosen from the list of single particle orbitals. A new pair of orbitals 
was formed via a simple unitary transformation 



with 6 sampled randomly from to 2ir. These moves were attempted 40 times each during 
a Monte Carlo pass (1000 times during the first pass). 

B. Addition/removal of kinks 

In this preliminary work only single and double "excitations" from the ground state were 
considered as candidates for kinks. That is, the difference between the ground state and 




(43) 



3=0 



4>'j = cos 9 <f)j + sin 9 & 
0^ = — sin 9 <pj + cos 9 <^ 



(44) 
(45) 
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an allowed excited state is the transfer of one or two electrons from occupied orbitals to 
unoccupied orbitals. The following scheme was therefore used for addition and removal of 
kinks. After sampling the rotations during the first Monte Carlo pass, the ground state was 
identified (this state did not change during the remainder of the simulation). A list of states 
corresponding to double and single excitations was constructed and used for the remainder 
of the simulation (the list would have been updated if the ground state had changed). A 
value of N g = 10 was used and at each Monte Carlo pass included an attempt to increase 
or decrease by 1 the value j in Eqn. EIBI 

C. Moving atoms 

Each Monte Carlo pass included an attempt to move each atom in turn, as in a standard 
Monte Carlo simulation. When an atom move was attempted, the single particle states were 
no longer orthogonal; the orbitals were orthonormalized using the Gram-Schmidt method. 
The step size for each move was 0.03 a.u. 

IV. SIMULATION DETAILS 

Monte Carlo simulations were performed using the sampling procedure described above. 
A temperature of l/fc^T = 3000 a.u. (~ 100K) was used; this was high enough to allow the 
relatively large geometric changes required to find the global minimum geometries, but not 
so high as to introduce large vibrational motion. All molecules were started in planar (and 
linear in the case of H 2 0) geometries to test the ability of the method to find the correct 
geometry. P = 3 x 10 10 and 1000 Monte Carlo passes were performed and averages were 
computed using the last 500 passes. One and two-electron integrals were evaluated using 
the C version programs included in PyQuante [25] and SCALAPACK[26] routines were 
used to perform matrix inversions. The 6-31G basis set was used and two simulations were 
performed for each molecule. In the first, no kinks were allowed resulting in a simulation 
using Eqn. A second simulation was performed using Eqn. H31 providing a simulation 
that included correlation. Calculations used 16 processors on the SuperHelix computer at 
LSU (www.cct.lsu.edu.) Correlation lengths of the energy and bond lengths ranged from 
50-150 Monte Carlo passes, reasonable values given the Monte Carlo step size. 
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V. RESULTS AND DISCUSSION 



The first molecule studied was H2O. Started in a linear geometry, the molecule quickly 
became bent and adopted the expected geometry. Fig. [T] shows the variation of total energy 
and ground state (HF) energy during the simulation using Eqn. H21 Fig- 121 displays the 
evolution of the different internuclear distances during the simulation. Several important 
features are evident from these figures. First, the electronic and geometric degrees of freedom 
evolve to their equilibrium values in a similar number of Monte Carlo passes. Second, 
the fluctuations in the Hartree-Fock energy are quite large (~.03 a.u. ~ 19 kcal/mol), 
demonstrating that the Monte Carlo procedure does not require a particularly accurate 
estimate for the Hartree-Fock ground state. Third, despite the fluctuations in the Hartree- 
Fock energy, the correlated energy has small fluctuations, which is to be expected of a 
good algorithm. From Fig. ^ and Table H] we can compare the correlated and uncorrelated 
methods used in this study. The Hartree-Fock estimator (Eqn. IT3*j) results in fluctuations in 
the energy estimator that are small and comparable to the fluctuations in the total energy 
estimator using Eqn. 1^31 

Tables lll lllll summarize the energies and geometries for all simulations. For comparison, 
we have calculated the K energies (with and without zero point energy correction) and 
geometries using Gaussian 98 (27 1 . 

Comparing the energies calculated using the path integral simulations with those obtained 
using Gaussian 98, we find that the path integral results using Qled are in close agreement 
with MP2 level energies and in much better agreement with MP2 energies than are the DFT 
results. The path integral energies lie above the K ab initio energies, as expected due to 
the entropy associated with the unitary transformations and finite temperature (classical) 
vibrational effects. The energies are below the zero point corrected ab initio energies due 
to the large zero point energy corrections in these molecules. The average geometries are in 
good agreement with ab initio results both in absolute bond lengths and angles and in the 
differences between Hartree-Fock and correlated (MP2/DFT) results, particularly in light 
of the ab initio bond lengths and angles being appropriate at K and the path integral 
results appropriate at 100 K. The largest difference between the path integral and ab initio 
geometries occurs in NH 3 , where the bond angles are larger by 2 - 3 degrees in the path 
integral calculations. Since the ab initio geometries are obtained at K and there is a 
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relatively low frequency vibrational mode in NH3, we performed a path integral calculation 
at a lower temperature to see if the average angles from a simulation came into better 
agreement with the K results. A simulation at (3 — 10000 a.u. (~ 30 K) resulted in 
energies and geometries shown in Table [HJ The bond angles are found to be in much better 
agreement with the K values. In addition, the energies are in better agreement with ab 
initio results. 

The Monte Carlo procedure selects excited states that make a significant contribution 
to the partition function. The number of possible 1- and 2- electron excited states range 
from 2240 for H2O to 5040 for CH4. The average number of excited states in the simulation 
ranges from 400-500. Thus, the Monte Carlo procedure is able to restrict the computational 
effort and bodes well for the scaling in larger systems. 

VI. CONCLUSION 

The kink-based path integral formulation has been extended to molecular systems. An 
approximate infinite order summation is used to include Hartree-Fock-like excited states in 
the ground state, correlated wavefunction. This procedure is necessary because all matrix 
elements are referenced to atom based primitive orbitals, which makes storage of the full N- 
electron density matrix too time consuming to be feasible. The estimator developed using 
this approach was compared to a Hartree-Fock-like method. In terms of geometries, the 
correlated method compares well with standard ab initio MP2 results (and are significantly 
better than DFT level results) and the Hartree-Fock-like geometries are in good agreement 
with temperature ab initio Hartree-Fock calculations. The treatment of geometric and 
electronic degrees of freedom on the same footing is a strength of this method. These initial 
results suggest this approach, combined with parallel computing, will provide an important 
alternative to standard ab initio methods, as well as the very successful Car-Parrinello DFT 
method. 

A direct comparison between the computational effort of conventional ab initio ap- 
proaches and the kink-based method is somewhat difficult since the bottleneck in a con- 
ventional simulation is the solution of the Hartree-Fock problem while in the path integral 
calculation a matrix inversion is the time consuming part of the calculation. It is possible, 
though, to discuss possible improvements to the kink-based approach. The time to invert a 
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matrix scales with the third power of the number of excited states, which (even with parallel 
computing) may prove to be a bottleneck to computational efficiency. However importance 
sampling using an importance function such as 



QT ED = J dK N j(* + YZffpT) (46) 
= Tr{W -(I + M )} 1 (47) 

rr = 4r^T ( 48 ) 



dxo 

r [ ( Xr .-l- F ) P I r apP 

Qled = dR N { i ° +1 -£; p (*o + r ^) P (49) 

reduces the computational effort significantly because only matrix multiplications are in- 
volved in the actual Monte Carlo moves. Some initial studies with this importance function 
showed a significant decrease in computational effort with only a minor decrease in precision. 

Also of interest is the question of size consistency. If a system is duplicated n times 
into n non-interacting systems, the partition function becomes the product of the individ- 
ual partition functions, which will guarantee size consistency. Also of interest from a size 
consistency point of view is whether the Monte Carlo estimator reaches this factorization 
limit. An examination of the leading terms in Eqn. EH indicates that r scales with n and 
r\ is independent of n. The latter feature can be understood in the following way. The 
denominators in Tq do not scale with n because the allowed excited states j are localized on 
one of the n systems. However, the derivative necessary to obtain I\ scales with 1/n. Since 
the number of excited states also scales with n, I\ will be independent of n. Therefore, the 
Monte Carlo estimator for n systems becomes 



x o + t^tJ w x ° exp ( ffrt 1 (50) 



This clearly is size consistent. 
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TABLE I: Energies, average number of excited states included in the path integral calculation 
(< N s >), and structural parameters for H2O. All energies and distances are in atomic units and 
numbers in parenthesis represent 2 standard deviations (95% confidence limits). <E> is the energy, 
including correlation, <Ehf > is the energy of the lowest energy state, < dun > is the average 
H-H bond length, < don > is the average O-H bond length, and < auou > is the average H-O-H 
angle. Path integral calculations were performed at (3 = 3000 a.u. ( ~ T = 100 K). Ab initio 



results were obtained using Gaussian 98 
correction (zpe). 



2jfl and are given with and without the zero point energy 



H 2 


<E> 


<Vhf > 


< N s > 


< dHH > 


< d H > 


< OLHOH > 


PI, Qhf (Eqn. CHI) 




-75.979(1) 





1.57(1) 


0.951(4) 


111(2) 


PI, Q LED (Eqn.El 


-76.096(2) 


-75.93(1) 


578 


1.57(1) 


0.968(4) 


109(1) 


ab initio HF (with zpe) 

(without zpe) 




-75.963 
-75.985 




1.57 


0.95 


112 


ab initio DFT (B3LYP, with zpe) 
(without zpe) 


-76.366 
-76.386 






1.58 


0.98 


108 


ab initio MP2 (with zpe) 

(without zpe) 


-76.092 
-76.113 






1.59 


0.97 


109 
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TABLE II: Energies, average number of excited states included in the path integral calculation 
(< N s >), and structural parameters for NH3. All energies and distances are in atomic units 
and numbers in parenthesis represent 2 standard deviations (95% confidence limits). <E> is the 
energy, including correlation, <Ehf > is the energy of the lowest energy state, < dun > is the 
average H-H bond length, < <Inh > is the average N-H bond length, and < cxhnh > is the average 
H-N-H angle. Path integral calculations were performed at (3 = 3000 a.u. ( ~ T = 100 K) except 
as noted. Ab initio results were obtained using Gaussian 98 27] and are given with and without 
the zero point energy correction (zpe). 



NH 3 


<E> 


<Vhf > 


<N S > 


< d H H > 


< dNH > 


< OtHNH > 


PI, Qhf (Eqn.CSJ 




-56.156(1) 





1.71(1) 


0.989(4) 


119(1) 


PI, Qled (Eqn.03) 


-56.240(2) 


-56.09(1) 


417 


1.70(1) 


1.000(3) 


117(1) 


PI, 0= 10000 a.u., Q L ed (Eqn.iSJ) 


-56.2760(2) 


-56.140(2) 


951 


1.694(2) 


1.009(1) 


114.2(2) 


ab initio HF (with zpe) 

(without zpe) 




-56.129 
-56.166 




1.68 


0.99 


116 


ab initio DFT (B3LYP, with zpe) 
(without zpe) 


-56.498 
-56.532 






1.71 


1.01 


116 


ab initio MP2 (with zpe) 

(without zpe) 


-56.244 
-56.280 






1.70 


1.01 


114 
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TABLE III: Energies, average number of excited states included in the path integral calculation 
(< N s >), and structural parameters for CH4. All energies and distances are in atomic units 
and numbers in parenthesis represent 2 standard deviations (95% confidence limits). <E> is the 
energy, including correlation, <Ehf > is the energy of the lowest energy state, < dun > is the 
average H-H bond length, < den > is the average C-H bond length, and < ancn > is the average 
H-C-H angle. Path integral calculations were performed at (3 = 3000 a.u. ( ~ T = 100 K) Ab 
initio results were obtained using Gaussian 98 2jfl and are given with and without the zero point 
energy correction (zpe). 



CH4 


<E> 


<Vhf > 


< N a > 


< dHH > 


< den > 


< OLHCH > 


PI, Qhf (Eqn. 03J 




-40.168(1) 





1.766(5) 


1.082(4) 


109.4(1) 


PI, Q LED (Eqn. Hi 


-40.254(4) 


-40.11(1) 


525 


1.782(6) 


1.092(4) 


109.4(4) 


ab initio HF (with zpe) 

(without zpe) 




-40.132 
-40.181 




1.78 


1.08 


110 


ab initio DFT (B3LYP, with zpe) 
(without zpe) 


-40.465 
-40.511 






1.79 


1.09 


109 


ab initio MP2 (with zpe) 

(without zpe) 


-40.233 
-40.279 






1.79 


1.10 


109 
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Figure Captions 



Figure 1. Energies during different Monte Carlo simulations of H2O. Energy during Qled 
simulation is the energy during a simulation using Eqn. 1321 Hartree-Fock energy during 
Qled simulation is the energy of the lowest energy state during a simulation using Eqn. H3J 
and energy during Qhf simulation is the energy during a simulation using Eqn. fTHl 
Figure 2. Internuclear distances during different Monte Carlo simulation of H 2 using 
Qled- The two hydrogen atoms are labeled Hi and H2. The initial values of the interatomic 
distances correspond to the initial linear geometry. 
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FIG. 2: 



23 



